#include <iostream>


extern "C" void computeGaussianGold(float delta, int radius);
extern "C" void bilateralFilterDepthGold(double* pSrc,
        double* pDest,
        float e_d,
        int w, int h, int r);
//variables
float gaussian[50];

void computeGaussianGold(float delta, int radius)
{
    for(int i = 0; i < 2 * radius + 1; i++)
    {
        int x = i - radius;
        gaussian[i] = exp(-(x * x) /
                          (2 * delta * delta));
    }
}

void bilateralFilterDepthGold(double* pSrc,
                              double* pDest,
                              float e_d,
                              int w, int h, int r)
{
    float domainDist, colorDist, factor;

    for(int y = 0; y < h; y++)
    {
        for(int x = 0; x < w; x++)
        {
            float t = 0.0f;
            float sum = 0.0f;
            float div = 1.0f /(2.0f*e_d*e_d);
            for(int i = -r; i <= r; i++)
            {
                int neighborY = y + i;

                if(neighborY> 0 && neighborY<h)
                {
                    for(int j = -r; j <= r; j++)
                    {
                        domainDist = gaussian[r + i] * gaussian[r + j];
                        int neighborX = x + j;
                        if(neighborY> 0 && neighborY<h)
                        {

                        float diff = fpSrc[y * w + x]-pSrc[neighborY * w + neighborX]);
                        colorDist = exp(-diff*diff*div);
                        factor = domainDist * colorDist;

                        sum += factor;
                        t += factor*float(pSrc[neighborY * w + neighborX]);

                        }
                    }
                }
            }
            pDest[y * w + x] = (double) (t / sum);
        }
    }
}
